Simulations of Blood Flow in Plain Cylindrical and Constricted Vessels with Single 

Cell Resolution 
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Understanding the physics of blood is challenging due to its nature as a suspension of soft particles 
and the fact that typical problems involve different scales. This is valid also for numerical investi- 
gations. In fact, many computational studies either neglect the existence of discrete cells or resolve 
relatively few cells very accurately. The authors recently developed a simple and highly efficient 
yet still particulate model with the aim to bridge the gap between currently applied methods. The 
present work focuses on its applicability to confined flows in vessels of diameters up to ~ 100 /im. 
For hematocrit values below 30%, a dependence of the apparent viscosity on the vessel diameter in 
agreement with experimental literature data is found. 
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I. INTRODUCTION 

Human blood can be approximated as a suspension 
of deformable red blood cells (RBCs, also called ery- 
throcytes) in a Newtonian liquid, the blood plasma. 
The other constituents like leukocytes and thrombocytes 
(blood platelets) can be neglected due to their low vol- 
ume concentrations. ^ Typical concentrations for RBCs 
are 40 to 50% under physiological conditions. In the ab- 
sence of external stresses, erythrocytes assume the shape 
of biconcave discs of approximately 8 /im diameter. S An 
understanding of their effect on the rheology and the clot- 
ting behavior of blood is necessary for the study of patho- 
logical deviations in the body and the design of microflu- 
idic devices for improved blood analysis. 

Well-established methods for the computer simulation 
of blood flow consist of an elaborate model of deformable 
cell membranes coupled to the surrounding plasma de- 
scribed by particle-based hydrodynamic methods , the 
lattice Boltzmann (LB) method 1^^, or the boundary in- 
tegral method 13. Others restrict themselves to a con- 
tinuous description at larger scales. 12) Our motivation 
is to bridge the gap between both classes of models by 
an intermediate approach that was published recently. ^ 
During the last decade, other groups already presented 
coarse-grained yet particulate models for blood at the 
mesoscale. Both Fenech et al. ti2J and Chesnutt and Mar- 
shall tiiJ developed discrete element models for the aggre- 
gation of large numbers of RBCs in two and three dimen- 
sions. These works, however, do not explicitly resolve 
the hydrodynamics of the blood plasma. Approaches 
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based on dissipative particle dynamics (DPD) directly 
model hydrodynamics. Boryczko et al. ^ applied their 
model that comprises plasma, deformable red cells, and 
capillary walls to the study of fibrin aggregation in the 
microcirculation. Filipovic et al. ^ employed DPD for 
the study of platelet adhesion to vessel walls and Pivkin 
et al. for the investigation of the effect of RBCs on 
the aggregation of blood platelets. Another simulation 
method that is well suited to account for hydrodynamics 
in confined geometries is the LB method which conse- 
quently has been applied in blood models of various levels 
of detail. Additionally to the models mentioned before 
already which are either based on a continuous fluid ^ or 
a suspension of fully deformable cells , there exist also 
more coarse grained particulate models employing the LB 
method: Sun et al. ^ model both red and white blood 
cells in two dimensions as rigid ellipses and spheres, the 
latter of which are, however, equipped with an elaborate 
wall-adhesion model. In the work by Hyakutake et al. l±£J 
the behavior of two-dimensional rigid spheres represent- 
ing RBCs at microvascular bifurcations is studied. 

Our coarse-grained blood model aims at a minimal res- 
olution of red blood cells which allows for a simple and 
highly efficient but still particulate description of blood 
as a suspension. Still, we do not give up explicitly mod- 
eling the suspending blood plasma or accounting for the 
non-spherical shape of RBCs. The ultimate goal is to per- 
form large-scale simulations that allow to study the flow 
in realistic geometries but also to link bulk properties, for 
example the effective viscosity, to phenomena at the level 
of single erythrocytes. Only a computationally efficient 
description allows the reliable accumulation of statistical 
properties in time-dependent flows which is necessary for 
this task. The improved understanding of the dynamic 
behavior of blood might be used for the optimization of 
macroscopic simulation methods. 
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The main idea of our model is to distinguish between 
the long-range hydrodynamic coupling of cells and the 
short-range interactions that are related to the complex 
mechanics, electrostatics, and the chemistry of the mem- 
branes. The short-range behavior of RBCs is described 
on a phenomenological level by means of anisotropic 
model potentials. I2J Long-range hydrodynamic interac- 
tions are accounted for by means of an LB method. liD 
Our model is well suited for the implementation of com- 
plex boundary conditions and an efficient parallelization 
on parallel supercomputers. Both are necessary for the 
study of realistic systems like branching vessels and the 
accumulation of statistically relevant data in bulk flow 
situations. In Section |n] and IIII| we briefly explain our 
approach. For an extended description we refer to our 
earlier publication. I2J Section IIVI contains — next to pa- 
rameter studies and a review of earlier results^ — an in- 
vestigation of the apparent blood viscosity in tubes and 
the related Fahraeus-Lindqvist effect. It closes with a 
qualitative view on constrictions and branching points in 
capillaries and data regarding the scalability of the code. 
In Section [Vj the conclusions from our work are drawn. 

II. HYDRODYNAMIC PART OF THE MODEL 

We apply a Bhatnagar-Gross-Krook (BGK) LB 
method for modeling the blood plasma. liSJ See for ex- 
ample the book of Succi for a comprehensive introduc- 
tion. liZJ The single particle distribution function n r (x, t) 
resembles the fluid traveling with one of r = 1, . . . , 19 
discrete velocities c r at the three-dimensional lattice po- 
sition x and discrete time t. Its evolution in time is de- 
termined by the lattice Boltzmann equation 

n r (x + c r , t + 1) = n r (x,i) — O (1) 

with 

n = n r (x,t) -<i(g(x,f),u(x,t)) ^ 

T 

being the BGK-collision term with a single relaxation 
time t. The equilibrium distribution function n° q (g,u) 
is an expansion of the Maxwell— Boltzmann distribution. 
e(x,i) = ^ r n r (x,t) and g(x, i)u(x, i) = £) r n r (x, t)c r 
can be identified as density and momentum. In the limit 
of small velocities and lattice spacings the Navier-Stokes 
equations are recovered with a kinematic viscosity of v = 
(2r — l)/6, where r = 1 in this study. 

For a coarse-grained description of the hydrodynamic 
interaction of cells and blood plasma, a method similar 
to the one by Aidun et al. modeling rigid particles of 
finite size is applied. [ 19 i 2 °l Starting point is the mid-link 
bounce-back boundary condition: the confining geometry 
is discretized on the lattice and all internal nodes are 
turned into fluid-less wall nodes. If x is such a node the 
updated distribution at x + c r is determined as 



with 

<(x,i) =n r (x,t)-n. (4) 

This corresponds to replacing the local distribution in 
direction r with the post-collision distribution n* (x, t) of 
the opposite direction f. To model boundaries moving 
with velocity v, Equation Q needs to be modified. The 
resulting update rule 

n r (x + c r ,t + 1) = n£(x + c r ,t) + C , (5) 

with 

2a c 

C = — ^£>(x + c r , t) c r v (6) 

was chosen consistently with n° q (g, u) for the general case 
u = v ^ 0. The lattice weights a Cr and the speed of 
sound c s are constants for a given set of discrete veloci- 
ties. The momentum 

A Pfp = (2n f + C)c f , (7) 

which is transferred during each time step by each sin- 
gle bounce-back process is used to calculate the result- 
ing force on the boundary. A freely moving particle i 
is modeled by such moving walls and is defined by its 
continuous position r^: all lattice nodes within a given 
theoretical particle volume — e. g. a sphere — centered at 
Ti are considered a wall and Equation is applied to 
links c r connected to its surface. When changes, even- 
tually the particle's discretization on the lattice needs to 
be updated. During this process, fluid nodes are created 
or vanish and the related change in total fluid momen- 
tum is balanced by an additional force on the respective 
particle. Still, the discretized particle shape undergoes 
fluctuations depending on the offset s of r 4 to the nearest 
lattice node. Their effect on hydrodynamics is quantified 
by measuring the translational and rotational drag coeffi- 
cient £ t and £ r of a sphere of radius R at a selection of dif- 
ferent offset vectors s. Since the model can be calibrated 
by assuming an effective hydrodynamic radius R* differ- 
ent from R , the drag coefficients are not normalized 
to the theoretical values but to their respective averages. 
The results in Figures Q] and [2] show that discretization 
effects depend on R but even for R = 2.5 typical de- 
viations are not larger than 10%. Typically, deviations 
in the rotational drag coefficient appear to be two times 
larger than the corresponding translational values. 

In contrast to the biconcave equilibrium shape of phys- 
iological RBCs and the previous tests with spherical par- 
ticles we choose a simplified ellipsoidal geometry that is 
defined by two distinct half-axes i?y and R± parallel and 
perpendicular to the unit vector 6^ which points along the 
direction of the axis of rotational symmetry of each par- 
ticle i. Since the cell-fluid interaction volumes are rigid 
we need to allow them to overlap in order to account for 
the deformability of real erythrocytes. We assume a pair 
of mutual forces 



n r (x + c r , t + 1) = n?(x + c r , t) , 



(3) 



F+=2<i(e,u = u)c r 



(8) 
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FIG. 1. Deviations of the translational drag coefficients £t 
from their average value £t at selected offsets s of the particle 
center to the nearest lattice node in the case of spherical 
particles with radius R. 



III. MODEL POTENTIALS 

In order to account for the complex behavior of real 
RBCs at small distances we add phenomenological pair 
potentials between cells. The idea is to develop simple 
model potentials and adjust their free parameters in or- 
der to match the pair interactions of cells. In a very 
similar way, the well-known Lennard- Jones potential is 
applied in classical molecular dynamics simulations to 
model atomic interactions. A fit can be achieved by com- 
parison of simulation results and experimental data from 
the literature, specially regarding blood rheology. For the 
moment, we concentrate on high shear rates 7 > 10s _1 
where aggregative effects play no role. S3 Therefore the 
model potential has to account only for deformation ef- 
fects. As a simple way to describe elastic deformability, 
we use the repulsive branch of a Hookean spring potential 
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FIG. 2. Deviations of the rotational drag coefficients £ r from 
their average value £ r at selected offsets s of the particle center 
Ti to the nearest lattice node in the case of spherical particles 
with radius R. 



and 



for the scalar displacement r,j of two cells i and j. With 
respect to the disc-shape of RBCs, we follow the approach 
of Berne and Pechukas^ and choose the energy and 
range parameters 
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(12) 

as functions of the orientations 6^ and 6j of the cells 
and their normalized center displacement f ij . We achieve 
an anisotropic potential with a zero-energy surface that 
is approximately that of ellipsoidal discs. Their half- 
axes parallel o\\ and perpendicular cr± to the symme- 
try axis enter Equation (fTTj) and (fT2f via a = 2a± and 
X = (<7jj — (i'j_)/((T^ + er^) whereas e determines the po- 
tential strength. For modeling the cell-wall interaction 
we assume a sphere with radius cr w = 1/2 at every lat- 
tice node on the surface of a vessel wall and implement 
similar forces as for the cell-cell interaction. Using 



pp 



2nf 1 (g, u = 0)c F 



pp 



(9) 



at each cell-cell link. This is exactly the momentum 
transfer during one time step due to the rigid-particle 
algorithm for a resting particle and an adjacent site with 
resting fluid at equilibrium and initial density g and thus 
compensates for the static fluid pressure. In case of close 
contact of cells with the confining geometry we proceed 
in a similar manner as for two cells but ignore forces on 
the system walls. 



c(6i,f ?x ) = 



1 - Xw (fixOi) 



(13) 



as a range parameter with a v 



and Xv 



(a 2 — CjJAo'ij + cr;^) allows to scale a potential with ra- 
dial symmetry to fit for the description of the interaction 
of a sphere and an ellipsoidal disc. ^ The parameter 
£(6^,6^) = e w remains constant in this case. fi X is the 
normalized center displacement of cell i and a wall node 
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FIG. 3. Two-dimensional cut as an outline of the model. 
Shown are two cells with their axes of rotational symmetry 
depicted by i>i/j- The volumes defined by the cell-cell inter- 
action are approximately ellipsoidal with half-axes Ox/il- The 
smaller ellipsoidal volumes (half-axes R±/\\) of the cell-plasma 
interaction are discretized on the underlying lattice. The cell- 
wall potential assumes spheres with radius <j w on all surface 
wall nodes x (taken from Reference^). 

Figure [3] shows a two-dimensional cut as an outline 
of the full model. For two RBCs the inner volume im- 
plementing the cell-plasma interaction with half axes i?y 
and R± is shown. Also depicted are the larger volumes 
that are defined by the range parameters au and <r± of 
the cell-cell and cell-wall interaction. A section of a vessel 
wall is visualized. The cell-wall interaction is based on 
the assumption of spheres with radius cr w at the location 
of the surface wall nodes. The forces and torques emerg- 
ing from the interaction of the cells with the fluid, other 
RBCs, and walls are integrated by a leap-frog scheme in 
order to evolve the system in time. Both the LB routines 
and the ones employed for treatment of the suspended 
cells use the same domain decomposition. For more de- 
tailed information concerning the model we refer to. 

IV. RESULTS 

All quantities can be converted from simulation units 
to physical units by multiplication with products of in- 
teger powers of the conversion factors Sx, St, and 5m for 
space, time, and mass. As a convention in this work, 
primed variables are used to distinguish quantities given 
in physical units from the same unprimed variable mea- 
sured in lattice units. Based on experimental measure- 
ments of RBC geometries, [2J the half-axes of the simpli- 
fied volume defining the cell-cell interaction in our model 
are set to 

a' ± = Afxm and a'« =4/3^m . (14) 

When choosing the spatial resolution quantified by the 
physical distance Sx corresponding to one lattice spac- 
ing, a compromise needs to be found: A low resolution 



severely reduces the computational effort but also reduces 
numerical accuracy. We decide for Sx = 2/3 /xm which 
allows a contiguous and closed volume of the cell-fluid in- 
teraction which still is completely comprised within the 
cell-cell interaction volume Equation (fT4"|) . On the ac- 
tual values of R± and R\\ we decide after studying their 
influence on the apparent viscosity /i a pp of the model sus- 
pension in plane Couette flow. 

Supposing that v matches the kinematic plasma vis- 
cosity of v' = 1.09 x 10~ 6 m 2 • s _1 , the time discretiza- 
tion is determined as St = 6.80 x 10~ 8 s. Since the 
largest shear rates employed in this work are of the or- 
der of 10 4 s _1 , sufficient temporal resolution is provided. 
Srn = 3.05 x 10~ 16 kg is chosen arbitrarily. All shear flow 
simulations reported here are performed on a system with 
a size of l x = 128 lattice units in x- and l y = l z = 40 lat- 
tice units in y- and z-direction or 85 x 27 2 /im 3 of real 
blood. Between the two yz-side planes a constant offset 
of the local fluid velocities in z-direction is imposed by an 
adaption of the Lees-Edwards shear boundary condition 
to the LB method. The apparent viscosity is calcu- 
lated from the average shear rate and shear stress. 12) Re- 
cently, an alternative method of viscosity measurement 
based on Kolmogorov flow was demonstrated which al- 
lows the employment of more simple periodic boundary 
conditions. ^ 

The influence of the volume of the cell-fluid interac- 
tion is studied as a function of the shear rate. Figure [4] 
shows the result for different Ru and R± at a cell num- 
ber concentration that varies by less than 5%. Generally, 
shear thinning is visible, but both the absolute viscosities 
and the change in viscosity per shear rate increase sig- 
nificantly for larger interaction volumes. Based on this 
and further parameter studies, 12J we choose with R' ± = 
11/3 //m the largest value investigated here and assuming 
fl||/o-|| = Rx/a± = 11/12 obtain R' = R'J'i = 11/9 /xm 
or Ru = 1.833 and R± = 5.5 measured in lattice units. A 
sphere with equal volume has a radius of R = 3.813. Fig- 
ures Q] and [5] suggest that at this resolution, we have to 
expect discretization errors of the order of some percents 
that are acceptable for our approach. 

In further simulations the effect of the stiffness param- 
eter of the cell-cell potential e is studied. The viscosity as 
a function of the shear rate increases with increasing e. 
For very stiff cells this dependence on the shear rate de- 
creases considerably which is in asymptotic consistency 
with the experimental results of Chien who measured 
the apparent shear viscosity of a suspension of artificially 
hardened RBCs and found a significantly increased yet 
mostly constant viscosity. At a cell-fluid volume concen- 
tration of 43% which seems sufficiently close to the hema- 
tocrit of 45% in the measurements done by Chien to 
justify a quantitative comparison, we find best agreement 
for e' = 1.47 x 10~ 15 J and use this parametrization for 
all following investigations. 12 

We now confine our model suspension in a cylindri- 
cal channel of diameter D and a length of 43 /im with 
periodic boundaries. Both cells and plasma are steadily 
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FIG. 4. Apparent viscosity ^ apP in dependence on the shear 
rate 7 as calculated for Couette flow for different sets of model 
parameters. a' ± = Aa'» = 4 /mi and e' = 1.47 x 1CP 14 J are 
kept fixed. While the volume implementing the cell-fluid in- 
teraction is varied from bottom to top as R'j_ — R'u = 1 /im, 
R' ± = 4R\\ = 8/3 fim, R' ± = 4R\ } = 10/3 /im, and R' ± = 
4Rm = 11/3 /im, the cell-fluid volume concentration <J> varies 

within 2.6% and 36%. Larger volumes 47ri?5_^||/3 lead to 
more pronounced shear thinning and generally higher viscosi- 
ties. 



driven by a volume force equivalent to a pressure gradi- 
ent dP/dz. We arbitrarily choose e' w = 1.47 x 10~ 16 J for 
the strength of the cell-wall interaction as a value that 
reliably prevents cells from penetrating the vessel wall. 
Figure \E\ visualizes the cells as the volumes defined by 
the cell-cell interaction and the channel wall as midplane 
between fluid and wall nodes in the case of D' = 63 /im 
and <f> = 25%. The pseudo-shear rate, which is defined 
via the volume flow rate Q, is v' = 4Q'/ (irD 13 ) = 61 s _1 . 
A preferential alignment of the cells largely perpendicular 
to the velocity gradient is visible. 

We compare radial velocity profiles for different flow 
velocities, a cell-fluid volume concentration of $ = 42%, 
and D' = 63 /im. While for high velocities the result 
looks parabolic in the central region, there is increas- 
ing blunting of the profile when the flow rate is reduced. 
The blunting can be understood as a consequence of the 
shear-thinning behavior of the model and is qualitatively 
consistent with experimental data from the literature. ^ 
Generally, apparent slip is visible close to the vessel wall. 
It is due to a cell depletion layer that to some extent can 
be controlled via e w . Thus, the observations can partly 
be described by an existing model assuming a homoge- 
neous core region with high hematocrit and consequently 
high viscosity and a cell-depleted boundary layer close to 
the vessel wall. I2i2§J 

As it is well known, the formation of this cell deple- 
tion layer influences the flow resistance of vessels which 
can be expressed in terms of their apparent viscosity. 
By inserting the measured volume flow rate Q through 
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FIG. 5. Steady flow through a cylindrical channel with a 
diameter of D — 63 /im at $ = 25% cell-fluid volume con- 
centration. The volumes defined by the cell-cell interaction 
are displayed at a cut parallel to the center axis. The flow 
is pointing from top to bottom and effects a pseudo-shear 
rate of v' = 61s _1 . Alignment of cells with the shear flow is 
observed. 



a cylindrical vessel into the theoretical expression for a 
Newtonian fluid 



Q = 



irD* dP 
128/zdz 



(15) 



and solving for fx, the respective apparent viscosity /i app 
can be determined. Except for Q, all known quantities in 
Equation (fT5f are constant and set as parameters of our 
simulation code. Only the flow rate undergoes stochastic 
fluctuations and — at the beginning of each simulation — 
shows a strong time dependence as the system relaxes 
from an arbitrary initial condition to a macroscopically 
steady state. Thus, typical simulations last for up to 
~ 10 7 time steps corresponding to ~ 1 s of physical time. 
Q is averaged over typically ~ 0.1s and its statistical 
error serves to estimate the accuracy of the resulting vis- 
cosity. The main dependency of /Lt app is on the hematocrit 
$. A comparison with experimental data in the case of 
Couette flow was published earlier. 12) However, /i app de- 
pends also on the vessel diameter D. This is known as 
Fahraeus-Lindqvist effect.^ Pries et al. combine a large 
set of experimental studies for v' > 50 s -1 and provide 
an empirical fit.l22J The resulting expression for /x a pp is 
plotted as a function of $ for four discrete values of D 
as lines in Figure [6] and [7] The symbols stand for simu- 
lation results at the same diameters D and pseudo-shear 
rates of (62 ± l)s~ 1 (Figure EJ) and (180 ± 5)s" 1 (Fig- 
ure [7]) . In consistency with the still significant shear- 
thinning that experiments but also our simulations Is 
exhibit above 7' = 50 s _1 , our results show a clear depen- 
dency on the pseudo-shear rate, which cannot be covered 
by the curves by Pries et al. that were obtained from aver- 
aging over different v 1 > 50s _1 .l22J Nevertheless, a com- 
parison confirms the agreement concerning the order of 
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FIG. 6. Dependence of the apparent viscosity /i apP in a cylin- 
drical vessel with diameter D on the volume concentration <E>. 
The lines show empirical results by Pries et al. 1221 for pseudo- 
shear rates v' > 50 s -1 with $ as tube hematocrit while the 
symbols represent our simulations for v = (62 ±1) s~\ Error 
bars are of the order of the symbol diameters and thus are not 
drawn. For $ < 0.3 there is good consistency of simulation 
results and experimental data. 

magnitude of our results and the observation that the ap- 
parent viscosity increases with D. This can be explained 
by the decreasing relative influence of the cell depletion 
layer. While the effect seems captured realistically for 
low volume concentrations $ < 0.3, the dependence of 
fi app on D becomes less clear for higher We explain 
this discrepancy by the fact that the thickness of the 
cell-depleted layer at high $ is determined by a balance 
of the short-range interactions of cells in the core act- 
ing towards a decrease of the depletion layer thickness 
and the short-range interactions of cells and the vessel 
wall acting towards its increase. Since the method pre- 
sented above models short-range interactions by means 
of potential forces and a constant pressure force Equa- 
tion © and velocity-dependent lubrication and lift 
forces can be covered only insufficiently in the case of 
very high volume concentrations. However, due to the 
Fahraeus effect, I2SJ <J> in smaller blood vessels is lower 
than the discharge hematocrit of typically 40-50%. Also 
the D-dependence found according to Pries et al. and 
displayed in Figure [Jj] and [7] is not very strong for the di- 
ameters studied. Thus the actual impact of the limitation 
of our model concerning <I> should not be overestimated 
when simulating vessels of comparable diameters. 

With a concluding study of the flow of nine RBCs 
through branching capillaries we demonstrate on a qual- 
itative level that our coarse-grained model is able to de- 
scribe problems where even single cells significantly af- 
fect the whole flow. While the diameter of the capillaries 
amounts to 9.3 fim, one of the branches features a steno- 



FIG. 7. Same plot as in Figure[6]but for v = (180±5) s _1 and 
two additional vessel diameters D. Due to shear-thinning, the 
viscosities here are slightly lower in general. For $ < 0.3, the 
dependence of /x app on D is captured again, but only on a 
more qualitative level. 



sis of only 5.3 /an. A cut of this setup is visualized in 
Figure [8] Both tube diameters and Reynolds numbers 
Re < 4 x 10~ 3 match physiological situations.^ We 
vary the cell-wall interaction stiffness e w and find that 
for e' w — 1.47 x 10~ 17 J, the erythrocytes easily pass the 
constriction as expected for healthy cells. tiJ The trajec- 
tories of the cell centers in this case are also shown in 
Figure [5] They visualize that — as expected from the lit- 
erature^ — a clear majority of cells is drawn into the 
unconstricted branch due to its higher flow rate. In con- 
sequence, the branch with the stenosis features a reduced 
hematocrit. Figure [9j which displays the time evolution 
of the relative volume flow rate in the constricted branch, 
makes clear that for e' w = 1.47 x 10~ 17 J this quantity 
changes only due to fluctuations induced by the mutable 
configuration of RBCs in the system. The constricted 
branch gets clogged for e' w = 1.47 x 10~ 15 J which is also 
reflected by Figure In the case of an intermediate value 
oi e' w = 1.47 x 10~ 16 J, cells are initially stopped but due 
to pressure fluctuations get squeezed through the steno- 
sis eventually. The acceleration by the cell-wall potential 
when leaving the constriction leads to the peaks of the 
flow rate visible in Figure [9j 

While the results for the straight channel studied first 
could still mostly be reproduced by a homogenous fluid 
with a specially tuned shear-rate or position dependent 
viscosity, this is not possible at the scale of capillar- 
ies. Here, our model allows to account for clearly par- 
ticulate effects like clogging or local changes of flow rate 
and pressure. In Figure [J2 it is demonstrated that such 
effects can lead to a distinct unsteadiness of the local 
cell volume concentration which is present also in human 
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FIG. 8. Branching capillaries visualized as a cut of the mid- 
plane between wall and fluid nodes. The channel diameter is 
5.3 ^tm at the stenosis in the upper branch and 9.3 /im other- 
wise. The flow direction is pointing from left to right. The 
plot also shows cell center trajectories at a cell-wall interac- 
tion stiffness of eL = 1.47 x 10" 17 J. 
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FIG. 9. Time evolution of the relative volume flow rate in 
the constricted branch of the system depicted in Figure [8] 
at different values of e' w . The clogging in the case of e' w = 
1.47 x 10~ 15 J becomes visible in a drop of the relative flow 
rate to less than 10%. While e' w = 1.47 x 1CT 17 J leads to a 
continuous flow situation (see Figure [8]), there are temporary 
drops and sharp peaks of the flow rate for e' w = 1.47 x 10 -16 J. 
They can be explained by RBCs being initially stopped at the 
stenosis and eventually squeezed through due to local pressure 
fluctuations (taken from Reference^). 



microvascular networks. I2£J 

Despite the simplifications of the model, parallel su- 
percomputers are necessary to simulate more realistic 
vessel networks or large bulk systems. This makes the 
scalability of the code crucial. For a quasi-homogenous 
chunk of suspension consisting of 1024 2 x 2048 lattice 
sites and 4.1 x 10 6 cells (see Figure [TP]) simulated on a 




FIG. 10. Schematic view of one of the square side planes of a 
benchmark system containing 1024 2 x 2048 lattice sites and 
more than 4 x 10 6 RBCs. The simulated volume resembles 
0.68 2 x 1.37 mm 3 of blood. 

BlueGene/P system, we achieve a parallel efficiency nor- 
malized to the case of 2048-fold parallelism of 85.6% on 
32768 and still 67.3% on 65536 cores. In comparison, the 
pure LB code without suspended RBCs shows a relative 
parallel efficiency of 98.0% on 65536 cores. The paral- 
lel performance of the combined code is mainly limited 
by the relation of the interaction range of a cell to the 
size of the computational domain dedicated to each task. 
The full strong scaling behavior for up to 65536 cores 
is shown in Figure 111! Compared to deformable par- 
ticle models, IS our method not only has a lower overall 
number of computations at a given resolution but is also 
easier to parallelize efficiently because each RBC has only 
six degrees of freedom. 

V. CONCLUSION 

We recently developed a computational model for the 
coarse-grained simulation of blood flow. I2J In the present 
work, additional bulk flow studies regarding the model 
parameters are shown. We further extend previous mea- 
surements of the apparent viscosity in cylindrical chan- 
nels to a comparison of different tube diameters between 
63 and 169 yum. While at higher hematocrit values the 
consistency remains limited to capturing the general be- 
havior of the viscosity and its order of magnitude, the 
Fahraeus-Lindqvist effect is reproduced very similarly to 
experimental data for $ < 0.3.122J With a qualitative 
study of unsteady flow and clogging in branching cap- 
illaries with a constriction, possible applications of our 
model are shown to range from bulk flow to capillary 
networks. Taking into account the good scalability of 
the code on massively parallel supercomputers, we be- 
lieve that our method can help to achieve a better un- 
derstanding of the links between the behavior of single 
cells and the properties of the encompassing larger flow 
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structures. 




2 4 8 16 32 
number of cores [1 024] 



FIG. 11. Strong scaling benchmark for the system depicted in 
Figure [TOI simulated on the BlueGene/P system at JSC. The 
relative speedup normalized to 1024 cores is plotted as a func- 
tion of the number of cores for the full model (plasma & cells) 
and a cell- free fluid volume of the same size (plasma). Due to 
higher memory requirements, the tasks in the combined case 
were performed by at least 2048 cores. 
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The simulation of hemodynamics is a very chal- 
lenging task. Recently, a scalable and simple particu- 
late model was developed that has the potential to bridge 
the gap that existing methods left open. Here, the model 
is applied to flows in vessels of diameters up to ~ 100 /im. 



